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ABSTRACT 

A new two-dimensional model of water flow in a 
hillslope has been implemented on the Massively 
Parallel Processor at GSFC. Flow in the soil both 
in the saturated and unsaturated zones, evaporation 
and overland flow are all modelled, and the 
rainfall rates are allowed to vary spatially. 
Previous models of this type had always been very 
limited computationally. This model takes less 
than a minute to model all the components of the 
hillslope water flow for a day. The model can now 
be used in sensitivity studies to specify which 
measurements should be taken and how accurate they 
should be to describe such flows for environmental 
studies . 


INTRODUCTION 

One important part of the global hydrological 
system is a catchment, which separates rainfall 
into evaporation, overland flow, and infiltration. 
For a heavy rain, infiltration excess reaches the 
stream first as overland flow. Part of the 
infiltrated water may then flow rapidly below the 
surface to re-emerge downslope or enter the stream. 
This is usually referred to as saturated subsurface 
flow. The rest reaches the unsaturated zone. The 
flow there is vertical and horizontal, and the 
latter component may eventually contribute to the 
stream flow. Another component which can 
contribute to the stream flow is horizontal flow in 
a perched water table above the bedrock. 

The primary output of catchment models is the 
hydrograph, in which the rainfall and fluxes to 
the stream from each of the above processes are 
plotted as a function of time. The rainfall 
rate and the sum of all the output fluxes are the 
usual data from a catchment, and a primary goal of 
catchment modelling is to understand the 
sensitivity of the output to the physical 
characteristics of the catchment, such as 
topography, cover type, soil characteristics, and 
antecedent moisture. 

Ref. 13 define catchment models as being of three 
basic types, but with overlapping characteristics 
so they may be considered a continuum. The first 
is stochastic. These models are statistical, in 
which time series of measured hydrographs (output) 
are correlated to rainfall (input) using classical 
time series analysis techniques. This leads quite 
naturally to parametric models, their second class, 
in which the parameters of the stochastic models 
are related empirically to the physical properties 


of the catchment. The third class contains 
deterministic models based on the laws of 
conservation of energy, mass, and momentum, usually 
expressed as time and space dependent differential 
equations. As these almost always contain non- 
measurable parameters which must be calibrated, 
deterministic models are partly parametric. 

There are many deterministic catchment models, but 
none of them includes all of the processes in the 
hydrological cycle. In part this is because we 
don't even know what they all are, due to the 
extreme complexity and variability of natural 
catchments. However, no existing model even 
includes all the processes previously 
described, because no serial computer can model 
them with a reasonable amount of computer time for 
a spatially variable catchment and for a long 
enough time period (Ref. 1,7,8,13.15). 

The concept of partial (or contributing) areas is 
one basis of our understanding of how catchments 
distribute rainfall (Ref. 17). Due to the spatial 
variability of catchment characteristics (soils, 
cover, topography), different areas handle the rain 
in different ways. For example, if the rain rate 
exceeds the infiltration capacity for a particular 
area, then the excess rain becomes overland flow. 
Once the soil is saturated, the water can flow 
rapidly below the surface and parallel to it. This 
process is referred to as saturated subsurface 
flow. The water will re-emerge somewhere 
downslope, adding to overland flow. The areas 
change over time, so the saturated partial area 
which contributes to overland flow varies in time 
as well as in space. 

We have tried to overcome the computing limitations 
by developing a model on the Massively Parallel 
Processor (MPP) . The model consists of a set of 
partial differential equations, solved in parallel, 
and so adapts naturally to a parallel architec- 
ture . The MPP hillslope model includes the 
following components: 

-- Surface retention 

- - A complete surface energy balance (tempera- 

ture and moisture) with separate evaporation 
rates from the soil, plants (with water 
extraction from the unsaturated zone) , and 
surface retention 

- - Overland flow 

-- Saturated subsurface flow parallel to the 
surface 
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-- Horizontal and vertical flow in the unsatur- 
ated zone 

Horizontal flow in an unconfined aquifer 

Our model is a vertical slice of a hillslope, so it 
is basically a two-dimensional model. It may be 
considered three-dimensional only if the gradients 
are all downslope, not across the slope. It is 
based on a catchment model of Ref. 11, which is 
simply a series of uncoupled one - dimens ional soil 
columns placed side by side. We have improved 
their design by allowing for horizontal flow in the 
unsaturated zone between the columns, and including 
the soil and surface temperatures. 

We decided at the beginning of this research effort 
to create one-, two-, and then three-dimensional 
models in succession. The one -dimensional model 
(Ref. 6) was compared to a similar one which runs 
on a serial machine (Ref. 5,10) to make sure the 
equations are solved correctly on the MPP, and as a 
timing benchmark. After the two-dimensional model 
is completely tested, we plan to develop a three- 
dimensional version. 

Our use of a parallel processor significantly 
reduces the execution time. Typically a 24 hour 
period may be modeled in about one CPU minute. 
Ref. 11 state that their model does not use 
excessive computer time on a serial machine, but 
they only present results from 6 hour simulations. 


THE TWO-DIMENSIONAL MODEL 


Saturated Zone 

The water table height in each column is HB. The 
horizontal flux is QB, and the vertical flux is QZ. 
The fluxes and vertical boundary conditions are 
calculated by the one dimensional Boussinesq 
equation (Ref. 14) . The flux into the water table 
from the unsaturated zone is modelled as the 
vertical hydraulic conductivity of the layer, and 
the bottom boundary condition is an input parameter 
representing an impervious layer or upward or 
downward seepage . 

The flux at the catchment divide is set to zero. 
At the seepage face the height HB is a fixed input 
parameter. Therefore the time derivation of HB is 
zero for the last column, and the discretized form 
of this derivative may be solved for the horizontal 
flux at the seepage face . This is the saturated 
zone flux which contributes to the hydrograph. 

Overland Flow 

If the surface water height is larger than a 
critical value, the overland flow flux is 
determined by Manning's equation (Ref. 7). 

The infiltration rate is basically the Green-Ampt 
model (Ref. 9), with the usual modification which 
replaces the depth of the wetting front with the 
cumulative infiltration: 

I(t) - a + b^/ t I(t') dt' ( 1 ) 

Surface Energy Balance 


The specifications for each of the components of 
the model given in the first section are described 
here as flux and continuity partial differential 
equations. The method of solution is also briefly 
described . 

Unsaturated Zone 

Moisture flow is modeled as described in Ref. 5, 
except we now have a horizontal component in the 
soil moisture flux. The surface temperature is 
modeled by the force - restore method. 

Boundary value fluxes must be specified for 
moisture at the top and bottom of the hillside 
(vertical direction) and at the hillslope divide 
and surfaces (horizontal direction) . The top 
boundary flux is the infiltration or evaporation 
rate, computed from the surface energy balance. 
The horizontal flux into the hillslope at the 
divide is zero. The horizontal flux at the 
hillslope surface depends on whether that cell is 
saturated. If it is and the sum of the vertical 
fluxes plus the horizontal flux into the cell from 
the interior of the hillslope would cause soil 
moisture to exceed saturation, then the flux onto 
the surface is set to whatever value is needed to 
keep moisture just as saturation. Otherwise, it is 
zero. This is the mechanism which allows 
subsurface return flow. 


The energy balance equation provides the surface 
fluxes : 

G-R + LE + LH ( 2 ) 

All fluxes are positive downward. G Is the heat 
absorbed by the soil, R is the net radiation flux, 
LE is the evapotranspiration energy flux, and H is 
the sensible heat. After finding the solution, the 
surface moisture flux q is set equal to the soil 
evaporation rate, and G is used in the force- 
restore model. The surface temperature needed to 
evaluate the fluxes is known from the force - restore 
equation. The latent and sensible heat fluxes are 
the usual resistance formulations. We imagine the 
soil and vegetation as one surface with the 
temperature T g . We also allow for some surface 
water storage. This affects the evaporation rates, 
because the surface resistance is zero for the 
fraction of the evaporation which comes from the 
stored water. 

Method of Solution 

The soil moisture and temperature continuity 
equations are solved by calculating the spatial 
derivatives of the moisture fluxes and then 
computing the time integral using numerical models. 

The soil is divided into cells by creating a grid 
of N layers and M columns of varying widths Az . and 
Ax i respectively, which are input parameters. i At a 
specified time the fluxes at the Interior 
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boundaries are calculated. The surface energy 
balance equations are evaluated and all boundary 
conditions applied. The continuity equations are 
of the form: 


dy -*■ 

f (t,y) 

dt 


(3) 


The vector y represents the state of the system in 
the unsaturated zone and ?(t,y) represents the 
model equations. This is solved with an Adams - 
Bashforth predictor -corrector method (Ref. 3,16). 
This solution is described in detail in Ref. 5. 
Since double precision is not available on the MPP, 
the form of the predictor-correction equations with 
the calculations done with the derivatives instead 
of the backward differences was used. New values 
of the state vector, y(p)(t+At) are predicted in 
terms of the previous derivatives. The derivatives 
are recalculated from the model equations, and then 
the corrected value of the state vector, 
y(c) (t+At) , is obtained. 

The difference between y(p) and y(c) is a reliable 
estimate of the discretization error, and the 
software determines if each element of this 
difference lies within a user- specif ied window. If 
all differences are smaller than this window, the 
integration step size (At) is doubled, leading to 
increased computational efficiency and reduced 
roundoff errors. If any difference is too large, 
the step size is halved. Doubling of the time step 
was accomplished by saving the previously 
calculated derivatives and using them. Thus, 
maximum accuracy could be retained. Where the time 
step could be doubled because the errors are small 
enough but there were insufficient back 
derivatives, doubling was postponed until there 
were sufficient back data. When the error window 
checks required that the time step be halved, three 
of the required derivatives for the predictor- 
corrector were available, and two were missing. 
The Runge-Kutta method was used to calculate these 
needed derivatives. The continuity equations for 
surface and saturated flow are solved using a 
Runge-Kutta method throughout. 


UTILIZING THE MPP ARCHITECTURE FOR SPEED 

Since identical calculations were needed at each 
soil cell, the mapping of the two dimensional model 
of the hi 1 Is lope was accomplished by assigning an 
individual processing element to each soil cell 
(see Fig. 1). Thus, the local memory of each 
processor contains the values which belong to that 
cell, i.e. moisture, position, thickness, depth, 
conductivities, etc. Surface temperature, deep 
soil temperature, cumulative infiltration, overland 
flow, and saturated flow were all stored as vectors 
in the same array as the moisture values since they 
were part of the state vector. 

The first step in the solution required calculation 
of the fluxes at the interior boundaries of the 



Figure 1. One processing element is assigned 
to one soil cell 

soil cells. These calculations involved only array 
arithmetic and nearest neighbor (in one direction 
for horizontal fluxes and in the other direction 
for vertical fluxes) calculations. Since the 
interconnect * scheme of the MPP is a nearest 
neighbor network, all of the array arithmetic and 
nearest neighbor calculations could be done in 
parallel. The next step in the solution required 
the surface energy balance equations be evaluated 
and the boundary conditions applied. These all 
involved vector calculations. Numerous input 
vectors were required to do these calculations over 
the course of a model run. Some were time 
dependent vectors such as the air temperature 
across the surface of the hillslope throughout the 
day and some were static throughout the model run, 
such as surface slope, surface roughness, and 
surface vegetation properties. These vectors were 
packed into array columns. To get the vector data 
to a convenient place to do calculations, the row 
and column broadcast capability of the MPP was 
used. This allows fast broadcast of one element 
from each row (column) to the other processor 
memories in the same row (column) (see Fig. 2). 

It is not necessary that the broadcast row (column) 
be composed only of elements in a horizontal 
(vertical) direction but merely that one element 
per column (row) be selected. The MPP's capability 
to select arbitrary areas of an array for 
calculation via boolean masks allowed the completed 
vector calculation results to be placed for example 
into the processor memories of only the surface of 
the hillslope. This combination of data movement 
via broadcast and boolean selection enabled the 
vector calculations to be done simply. In 
addition, since many of the vector calculations 
were similar, it was possible to do more than one 
set at a time. 

Once the derivatives were calculated, the 
predictor - corrector equations were used and the 
differences between them found. The tests on the 
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MODEL OUTPUT 


Figure 2. The row and column broadcasting 
feature of the MPP allows quick 
movement of data for vector 
calculations 

halving (doubling) converted to a hardware 
instruction on the MPP and could thus be done in 
parallel. This global testing ability of the MPP 
was also used to decide if whole blocks of code 
needed to be executed or could be skipped. This 
occurred for example with the infiltration 
calculations under surface saturation. If no part 
of the surface was saturated, then these 
calculations could be skipped entirely. This also 
contributed to the overall speed of execution. 

In summary, the program's speed was achieved 
through array arithmetic (masked and unmasked), 
parallel data movement through nearest neighbor 
communication and row and column broadcasting, and 
global testing of conditions using 'any' or 'all' 
for the purpose of choice in the next set of 
calculations. All of these fitted naturally with 
the MPP architecture and the computational 
requirements of the model. A comparison of the 
times (see Table 1) for the model as it has evolved 
from a 14 layer, one -dimensional limited flow model 
to the current two dimensional model shows that a 
single day of data run through the model requires 
only about a minute of CPU time . 

Table 1. Timing measurements comparing MPP and 

a serial processor for 24 hours of data 
processed . 

One Dimensional Model 

(14 soil layers, no rain, vertical flows only) 

IBM (Full processing capability): 4 sec 

MPP (14/16384 processors): 10 sec 

Two Dimensional Model 

(102 soil layers, 102 soil columns, horizontal 
and vertical unsaturated flows, saturated flow, 
overland flow, one hour of rain) 

MPP: 57 sec 


We have not yet completed unit testing of all the 
processes in the model. Here we present the 
results of one test, which includes the surface 
energy balance of and infiltration into an 
initially very dry sandy loam soil. 

The hillslope is divided into 102 columns of width 
.5 meters each. The first column has 100 soil 
layers of thickness .1 m and the bottom two layers 
.5 meters. The last column has only the bottom two 
layers. The slope is a line drawn from the top of 
the first column to the top of the last, so the 
area modeled is a right triangle with height 11 
meters and base 60 meters. These soil cells plus 
the additional cells for temperature, infiltration, 
overland flow and saturated flow use approximately 
one -third of the Array Unit Processor capacity. 

The initial volume t^ic_ jio is ture in the unsaturated 
zone is set to .05 m m everywhere. To model a 
sandy loam we have set the parameters in the 
hydraulic conductivity and matrix potential models 
to 0 s - .375, K s « 2.8 x 10 ms - -.43 m, 

and b - 5. These values were derived from fits to 
the characteristic curves measured during an 
experiment near Phoenix in 1972 (Ref. 12). They 
were reused for each of the 6 days modeled here. 
Ref. 4 show how these data were fitted to the 
surface energy balance model. The rainfall rate 
was 1.6 cm h for the first 3 hours. 

Perhaps the most important result is that the 
simulations took approximately 1 minute of CPU time 
per 24 hour period, or 6 minutes for the entire 6 
day run. In numerical simulations on earth science 
problems, computer runs of an hour or more are not 
uncommon. In such a time, it is feasible to 
simulate 2 months or more of model time on the MPP. 
This will allow for simulations of many storms and 
inter-storm periods. 

Figure 3a and 3b show the force restore solutions 
to the surface and deep soil temperatures as 
functions of time and column number. Time zero is 
the start of the simulation, which here Is 
midnight. Column 1 is at the hillslope divide 
and column 102 is at the seepage face. It is 
difficult from these plots to project the daily 
maximum value onto the time axis, but for each day 
this occurs at 2 p.m, The temperatures range from 
22 to 40 ( C) , increasing as the soil surface 

dries. The temperatures in the last three columns 
show some problems, which we are examining. 

Figure 4a shows soil moisture in the top soil layer 
as a function of time and position. The rapid rise 
as the initially dry soil absorbs all the rain and 
the subsequent decline over the next 5 days are 
physically realistic. 

Figure 4b shows the soil moisture profile in column 
50 (halfway down the hillslope) as a function of 
time. This shows that the moisture never 
penetrates deeper than the top 5 layers, or .5 
meters. It also shows that after 2 days the 
surface exhibits small oscillations about a value 
of .05 (same as in Fig. 4a), increases to a value 
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Figure 3a. 


Force - res tore solution for the surface 
temperature as a function of time and 
position on the hillslope 


Figure 4a. Surface soil moisture as a function 

of time and position on the hillslope 



Figure 3b. Force-restore solution for the deep 
soil temperature as a function of 
time and position on the hillslope 


of about .12 at about .3 meters then 
an unchanging value of .05 below .5 
the dynamic zone seems to be the top 


decreases to 
meters. Thus, 
. 5 meters . 


Figure 4c shows the variation of 
moisture as a function of time . 
infiltration and evaporation 
capillary action, can be seen. 


the top cell soil 
The effects of 
as well as of 


Figure 4b. Soil moisture profile for column 50 as 
a function of time 

Figure 5a shows the infiltration rate as a function 
of time and. position. The maximum rate shown here 
(4 4 x 10 cm s ) equals the rain rate, 1.6 cm 
h" 1 . Figure 5b shows the cumulative evaporation 
everywhere as 4.8 cm, exactly equal to the 
cumulative rainfall. For this simulation then 
all the rain immediately infiltrated into the soil 
surface Figure 5b also shows that the cumulative 
infiltration calculation is correct. There is no 
surface retention. 
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radiation, a rather large value. The problem is 
not in the values for thermal conductivity and heat 
capacity, as may be seen in Figure 7. These vary 
with soil moisture as they should. 



^4 



Figure 5a. Infiltration rate as a function of 
position and time 

The surface energy balance fluxes are plotted in 
Figs. 6a-6d. The net radiation (Fig. 6a) is the 
data used to drive the energy balance model. These 
are the same very day, as we simply reused the 24 
hour data set each day. The latent heat flux (Fig. 
6b) decreases each day as the soil dries out. The 
sensible heat flux (Fig. 6c) exhibits peculiar 
behavior, being predominantly positive (towards the 
soil in the sign convention of Eq . 2) for the first 
4 days and negative thereafter. Finally, Figure 6d 
shows the soil heat flux. It is positive during 
the day as it should be for a soil surface which is 
getting warmer every day, but it is also 50% of net 


Figure 5b. Cumulative infiltration as a function 
of position and time 



Figure 6a. Surface net radiation as a function of 
time and position 

These peculiarities in the surface fluxes are most 
likely due to the use of the same net radiation 
every day, which cannot be representative of all 
the surface conditions modeled here. This is being 
checked out by using modeled instead of measured 
radiation . 
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Figure 6c. Sensible heat flux as a function of 
time and position 


SUMMARY 

We have presented a new model of the hydrological 
response of a hillslope to rain. It runs on a SIMD 
parallel architecture computer, the Massively 
Parallel Processor, at Goddard Space Flight Center. 
Its major advantage over other models of its type 
is its much reduced execution times (due to the 
parallel architecture of the MPP) from what one 
gets on a serial machine. This allows the model to 
include more of the hydrological processes than any 
other model has been able to, including saturated 
subsurface flow and a sophisticated surface energy 
balance . 


Figure 6d. Soil heat flux as a function of 
time and position 


^3 



Figure 7a. Thermal conductivity of the top soil 
layer as a function of time and 
position 
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Figure 7b. Heat capacity of the top soil layer 
as a function of time and position 
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